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ABSTRACT 

Oscillatory double-diffusive convection (ODBC) (also known as semi¬ 
convection) refers to a type of double diffusive instability that occurs in regions 
of planetary and stellar interiors which have a destabilizing thermal stratihcation 
and a stabilizing mean molecular weight stratihcation. In this series of papers, 
we use an extensive suite of three-dimensional (3D) numerical simulations to 


quantify the transport of heat and chemical species by ODBC. Rosenblum et ah 


(2011) hrst showed that ODBC can either spontaneously form layers, which 


signihcantly enhance the transport of heat and chemical species compared to mi¬ 
croscopic transport, or remain in a state dominated by large scale gravity waves, 
in which there is a more modest enhancement of the turbulent transport rates. 
Subsequent studies in this series have focused on identifying under what condi¬ 
tions layers form (Mirouh et ah [20 12), and quantifying transport through layered 


systems (Wood et ah 2013). Here we proceed to characterize transport through 


systems that are unstable to the ODBC instability, but do not undergo spon¬ 
taneous layer formation. We measure the thermal and compositional fluxes in 
non-layered ODBC from both 2D and 3D numerical simulations and show that 
3D simulations are well approximated by similar simulations in a 2D domain. 
We hnd that the turbulent mixing rate in this regime is weak and can, to a hrst 
level approximation, be neglected. We conclude by summarizing the hndings of 
papers I through III into a single prescription for transport by ODBC. 


Subject headings: convection, hydrodynamics, planets and satellites: general, stars: 


interiors 
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Introduction 


1.1. Background 


Stratified fluids which have a stabilizing compositional gradient (Ledoux stable) and 
a destabilizing thermal gradient (Schwarzschild unstable) are common in a variety of 


astrophysical objects. This was hrst discussed by Schwarzschild & Harm (1958) in the 


context of the growing convection zones of massive stars (> lOM^^^). However, Walin 


(1964) and Kato (1966) were the hrst to realize that huid instabilities may develop in such 


an environment, driving the turbulent transport of heat and chemical species. They are 
now commonly referred to as oscillatory diffusive instabilities, and they saturate into a 
non-linear weakly turbulent regime called Oscillatory Double Diffusive convection (ODDC). 
This nomenclature derives from the oscillatory nature of the basic unstable huid motions 


which bear resemblance to over-stable internal gravity waves (Kato 1966). The efficiency 


of mixing (of heat and chemical species) in ODDC is potentially signihcant to evolution 
models for stars (Langer et al.|[l9^; Merryheld 1995) and giant planets (Stevenson 1982). 


The conditions for ODDC to occur are dehned by three non-dimensional parameters 


(Baines & Gill 1969). The Prandtl number, Pr, and the inverse Lewis number, r, are 


dehned as 


Pr = ^ 

KT 


q- = In. 

Kt 


( 1 ) 


where v is the viscosity, and kt and are the thermal and compositional dihusivities, 
respectively. The third parameter, the inverse density ratio, is dehned as 
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( 2 ) 


where T, p and p are the temperature, pressure and mean molecular weight, and where 


5=1^1 and(^=|^ 

^ din PL 
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are thermodynamic derivatives of the equation of state. The 
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subscript “ad” is used to denote an adiabatic gradient (a derivative at constant entropy). 


ODBC occurs when ^ is within the following range (iKato 19661 Walin 19641), 




( 3 ) 


Pr + r 

When is less than one, the system is unstable to overturning convection, and when R^^ 
is greater than R~^ the system is completely stable. A useful proxy for the inverse density 
ratio is the so-called “reduced inverse density ratio” parameter ( Mirouh et al.||2012 ), defined 
as, 

Rjl - 1 

( 4 ) 


r = 


Rq" - 1 

Pr+l _ 1 
Pr+T 


This parameter maps the entire ODBC range to the interval [0,1] where 0 marks the onset 
of overturning convection, and 1 marks the boundary between the OBBC parameter space 
and that of complete stability. 


1.2. Recent studies and current work 

In previous papers of this series we have set out to systematically characterize the 
different types of behavior demonstrated in OBBC. By analyzing 3B direct numerical 


simulations, Rosenblum et ah (2011) hrst identihed two general classes of behavior: one in 


which thermo-compositional layers emerge after the growth and saturation of the primary 
OBBC instability (and cause a significant increase in transport), and one in which layers 
do not form and where the dynamics are dominated by what they called “homogeneous 


turbulence”. Mirouh et ah (2012) built on the work of Rosenblum et ah (2011) by 


developing a semi-analytical rule for identifying the regions of parameter space where 
layers do and do not form. They also determined that non-layered OBBC is ultimately 
dominated by large-scale internal gravity waves which (surprisingly) also augment thermal 


and compositional transport, though not as much as in the layered case. Wood et ah (2013) 
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then studied the dynamics and transport properties of layered ODBC. In this work we now 
focus on characterizing the behavior of non-layered ODBC and quantifying the associated 
thermo-compositional mixing it induces. 

In Section we present our mathematical model for the dynamics of OBBC. We also 
briefly describe how we study it using BNS, and discuss the metrics by which we measure 
thermal and compositional fluxes. In Section we describe the evolution of a typical 
non-layered simulation and discuss the difficulties of running numerical simulations in 
extreme parameter regimes. In Section we discuss the effectiveness of 2B simulations for 
modeling OBBC compared to the full 3B BNS. We present the results of our numerical 
experiments in 2B and 3B, focussing in particular on the measurements of thermal and 
compositional fluxes in Section Finally, in Section we present our conclusions, and 
summarize the hndings of papers I, II, and III of this series. 


2. Mathematical model and numerical simulations 


In what follows, we consider a region that is signihcantly smaller than the density 
scale height of a typical star or planet, and where flow speeds are much smaller than the 


sound speed of the medium. This allows us to use the Boussinesq approximation (Spiegel & 


Veronis||1960 ) and to ignore the spherical geometry of an actual stellar or planetary interior. 
We consider a Cartesian domain where gravity is oriented in the negative 2 ;-direction, and 
we assume constant background gradients of temperature and mean molecular weight that 
are dehned as follows. 


_dT_Tdp 
dz p dz ' 
dfi fJ'dp 
dz p dz ^' 


( 5 ) 




We use the following linearized equation of state, 


P- = -af + PfL, (6) 

Po 

where p, T and fl are dimensional perturbations to the background profiles of density, 
temperature and composition, and where po is the mean density of the region. The 
parameters a and f3 are the coefficient of thermal expansion and coefficient of compositional 
contraction, respectively, and are defined as. 
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P 


(7) 


Using these assumptions, the standard non-dimensional governing equations for ODBC 


'Rosenblum et ah 

2011 

Mirouh et ah 

2012 


va 


1 

Pr 



u ■ Vu 


-h u ■ VT — w 
ot 


dp 


-|- u ■ V/i — Rq 


0 , 

-Vp + {f - p)ez + V^u, 

V^T, 

rV'h, 


( 8 ) 


where u is the velocity field, and all values with tildes are now non-dimensional. To arrive 
at these equations, we define the unit length in terms of (dimensional) constants. 


[/] = d 



(9) 


where g is the gravitational acceleration, and where the parameter Tq^ is the background 
adiabatic temperature gradient which is defined as. 
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Using these constants, we also construct the units of time, [t], temperature. 


T 


and 


mean molecular weight, [/i], as 


[t] 



[/i] 


K'j' 

d\To.-T^\. 


( 11 ) 


Furthermore, we can also rewrite the inverse density ratio Rq ^ in terms of the 
parameters introduced above as. 


R 


-1 

0 


Iho^ I 

“|roz -ifcl ’ 


( 12 ) 


As the physical characteristics of stellar and planetary interiors cannot be easily 
reproduced in laboratory experiments, we must rely on numerical simulations to gain insight 
about the nonlinear behavior of ODBC in these objects. We solve the set of equations in 


using a pseudo-spectral code (Traxler et ah 2011) where all perturbations are triply 
periodic in the domain. Much of the experimental data presented in subsequent sections 


was generated by Rosenblum et ah (2011) and Mirouh et al. (2012) in simulation runs using 
this code. We have generated new data specihcally for this work as well, including a series 
of 2D simulations which we will compare to the 3D ones. 


The main quantities of interest we extract from these numerical simulations are the 
turbulent vertical fluxes of temperature and chemical species, {wT) and (tc/i), respectively, 
where the angle brackets represent a spatial integral over the entire computational domain. 
It is often useful to express these fluxes in terms of thermal and compositional Nusselt 
numbers (Nut and Nu^) which are the ratios of the total fluxes to the diffusive fluxes (of 
temperature and composition). The turbulent fluxes are dehned as follows in terms of Nut 
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and Nu^, 





r^turb 
tj. — 

{wT) = Nut — 1, 



r^turb _ 

(wp) = rRp^ (Nu^ - 1) . 

(13) 


As shown by Wood et al. (2013), these fluxes often exhibit fast oscillations with large 


amplitudes due to the gravity waves, so for the purposes of analysis, it can be more useful 
to consider related quantities known as the thermal and compositional dissipations, {|VTp) 
and {|V/ip). Indeed, taking a spatial integral of the thermal and chemical evolution 
equations, and then assuming that the system is in a statistically stationary state, we And 
that the dissipations are related to the turbulent fluxes and Nusselt numbers by. 


NuT-l = {h;T) = (|VT|2), 
{wp) (|V/1|2) 


Nu,, 


1 = 




(14) 


where the bars indicate temporal averages over the entire statistically stationary period. 

In practice, these are good approximations even when the temporal averaging is done over 
short periods, so in what follows we assume similar relations between the instantaneous 
Nusselt numbers and dissipations as well. This is advantageous because the dissipations are 
less sensitive to the oscillations of gravity waves than the fluxes, and are therefore easier to 
analyze. 


3. Behaviors of ODDC 

In what follows, we show the results of a typical ODDC simulation. The simulation 
presented has Pr = r = 0.03 and = 7.87, and was run at an effective resolution of 
192^ (the simulation domain has dimensions (lOOd)^). We first describe these results purely 
qualitatively, then move on to a more quantitative analysis. 
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3.1. Qualitative study 


Figure [T] shows the vertical component of the velocity held at early times, and Figure]^ 
shows the y-component of the velocity held later on in the simulation. At very early times, 
we hrst see the development of the fastest growing modes of the basic ODBC instability. 


which resemble tubes of vertically oscillating huid (shown in Figure la). This primary 
growth phase ends when the basic instability saturates due to nonlinear interactions 
inherent to the problem (see Figures [lb|. 

After the primary saturation, the small-scale fastest growing modes of the primary 
instability stop growing, but other larger-scale modes slowly continue to gain energy 


(amounting to a secondary phase of growth). From Figure 2a we see that the latter (which 
ultimately come to dominate the energetics of the system) generally have the largest possible 
scale in the horizontal direction. As the system evolves, energy is funneled into modes of 
progressively larger vertical scale until the secondary growth phase saturates and the system 
reaches a statistically stationary state. The dynamics of this state are characterized by 
gravity wave oscillations on fast timescales whose amplitudes are modulated chaotically and 
intermittently. This intermittency appears to be caused by nonlinear interactions between 
large scale gravity modes and large scale horizontal shearing modes. Indeed, we regularly 
observe the emergence of a strong horizontal shear layer that temporarily suppresses the 
wave field (Figure [2bl). The shear then dissipates, and the system proceeds as before. Figure 
shows the distinct differences in the y velocity field between a gravity-wave-dominated 
phase and a shear-dominated phase. 
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3.2. Quantitative study 


In order to study this system in a more quantitative way we now investigate the energy 
contained in individual modes. We shall refer to specihc spatial modes by the number 
of wavelengths in the x, y, and z directions. So mode {l,m,k) would refer to a mode 
with horizontal wave numbers ^ and (in the x and y directions), and vertical wave 
number ^ (in the z direction). We can quantify the transfer of energy to larger scales by 
considering the amount of energy in a mode “family”. A family of modes consists of all 
the modes with equivalent spatial structures given the symmetries between the x and y 
directions in the domain, and negative and positive wave numbers in each spatial direction. 
For example, the mode family 102 contains the modes (1,0,2) , (-1,0,2) , (0,1,2) , (0,-1,2) , 


(l,0,-2) , (-l,0,-2) , (0,l,-2) , (0,-l,-2) (Traxler et ah 2011). 


Consistent with the qualitative results in Figure Figure shows that at early times, 
the dynamics of ODBC are dominated by horizontally small scale modes {y/W+^iv? 8), 

with no structure in the vertical direction (with vertical wavenumber k = 0). Around 
t = 2500, the primary instability saturates. [Mirouh et ah (2012) demonstrated that the 
level at which the primary instability saturates can be used to identify regions of parameter 
space where layer formation is expected to occur. However, the primary saturation level is 
of little use for characterizing the long term transport properties of our non-layered ODBC 
because a secondary growth phase occurs after primary saturation, augmenting the thermal 
and compositional fluxes. From Figure]^ we see that while the fastest growing modes of 
the primary instability stop growing at saturation, the mode family (1,0,5) continues to 
grow, and for a brief time becomes the most energetic mode family in the system. As time 
goes on, however, the mode family (1,0,4) supplants (1,0,5) as the most energetic, which is 
in turn overtaken by the mode family (1,0,3). For this particular simulation, the handoff 
of energy to larger scale modes stops with mode family (1,0,3); mode family (1,0,2) never 
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becomes dominant. Crucially, Figures and also reveal the growth of the energy in 
large scale shearing modes with purely horizontal fluid motions (mode families (0,0,1) and 
(0,0,2)). This is unexpected because these modes are not directly excited by the primary 
ODBC instability. Instead their growth must arise from nonlinear interactions between 
rapidly growing ODBC modes. 

Around t = 5000, after mode family (1,0,3) becomes dominant, the secondary growth 
phase appears to end, saturating into a statistically steady turbulent state. However, Figure 
1^ shows periodic bursts of growth in the shearing mode energies, suggesting intermittent 
(yet powerful) interactions between the shearing modes and the dominant gravity wave 
modes. In these interactions, illustrated in more detail in Figure the growth of gravity 
waves excites the rapid growth of horizontal shearing modes, which in turn causes a decay 
in the amplitude of the large scale gravity waves. Without the wave held to amplify it, the 
shear then decays viscously. This hnally allows the energy in the gravity waves to ramp 
back up again. While this interaction does not always manifest itself as such a distinct 
sawtooth oscillation, it is still present in one form or another in all gravity-wave-dominated 
OBBC simulations. Figure also shows that the interaction between the shear and gravity 
waves also leads to intermittent modulation of the thermal and compositional huxes. 

While the intermittency in the huxes caused by the shear is interesting in its own right, 
for the purpose of modeling OBBC transport in planetary or stellar evolution calculations, 
we are more concerned with estimating the mean huxes over signihcant periods of time. 
These mean huxes at secondary saturation depend on the parameters of the system {Rq^, 
Pr, and r). The results shown in this section, which were obtained at moderate and 
moderate Pr and r, suggest that turbulent transport in non-layered OBBC is weak. To 
conhrm this we need to run numerical experiments at larger i?g ^ and smaller Pr and r. 
Probing this region of parameter space is difficult, however, because 3B simulations at low 
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Pr and r can be compntationally very expensive, particnlarly for valnes of that are 
close to marginal stability [Rq^ —)■ p^)• The small valnes of Pr and r lead to small-scale 
tnrbnlent featnres with steep gradients of velocity, temperatnre, and composition, which 
necessitate high spatial resolntion. Fnrthermore, a larger Rq^ leads to higher freqnency 
oscillations of the basic ODBC modes, necessitating smaller time steps. Given these 
challenges, in the next section we discnss the possibility of nsing 2D ODBC simnlations as 
a potential snrrogate for fnll 3D simnlations at these extreme regions of parameter space. 


4. 2D simulations 


Simnlations of 2D ODBC systems are compntationally inexpensive and are also less 
intensive in terms of data storage than 3D simnlations. For this reason, we have rnn a 
series of tests to compare both qnalitative behavior and qnantitative estimates of the flnxes 
(and other system diagnostics) in 2D and 3D. Fortunately, as we see from Figure the 
secondary saturation level in the 2D simulation at Pr = r = 0.03 and R^^ = 7.87 is very 
similar to the 3D simulation analyzed in the previous section. This is, in fact, generally 
the case for each parameter regime where we have both 2D and 3D simulations, which is 
surprising in light of the fact that 2D simulations of overturning convection (i.e. Rq ■'<1) 


behave very differently from their 3D counterparts (Schmalzl et ah 2004; van der Poel et ah 


2013). In the results that follow in section]^ we shall rely heavily on 2D simulations to 


draw conclusions about turbulent fluxes through non-layered ODBC. 


Measurements of mode family energies show that key physical processes that dictate 
the behavior of 3D ODBC simulations are present in the 2D simulations as well. Figure 
[^explores the energetics of the gravity waves and shear, and shows that the fractions of 
energy in each type of mode are of the same in order in both cases. This is important 
because together these two types of modes contain most of the energy in non-layered 









systems after secondary saturation. 


The computational economy of 2D simulations makes other types of analysis easier as 
well, such as running simulations in larger domains. In the previous section, we showed that 
after primary saturation the dominant gravity wave modes have horizontal wavelengths 
commensurate with the domain size. Also we showed that energy is transferred to modes 
with progressively larger vertical wavelengths. This raises the question of whether this 
energy transfer would always terminate at a vertical wavelength that depends on the 
domain size. For example, will the dominant mode after secondary saturation in a (200(i)^ 
domain have a vertical wavelength that is twice that of the dominant mode in a (lOOd)^ 
domain? More importantly, do the fluxes depend on the domain size? 


Using 2D data we hnd that in all but one case, doubling the domain size in each 
direction, leaves the vertical wavelength of the dominant mode unchanged. By contrast, 
the horizontal wavelength of the dominant mode always grows to the largest possible scale 
allowed by the domain. As a consequence, the dominant modes in the larger boxes are 
inclined more toward the horizontal than in the smaller ones. Importantly though. Figures 


6a and 6b show that despite some qualitative differences between simulations with domains 


of different dimensions, we hnd that the time-averaged huxes of temperature and chemical 
composition do not depend strongly on box size (they are within ~ 10% of one another). 


5. Results and Discussion 

We now analyze the results of all numerical experiments done in the 2D and 3D 
computational domains. We evaluate thermal and compositional huxes in terms of the 
Nusselt numbers, which we calculate from thermal and compositional dissipation data, as 


described in (14). More specihcally, the quantities we consider are Nu^ — 1 and Nu^ — 1 
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which, conveniently, are also measnres of the tnrbnlent flnxes in nnits of the diffusive fluxes. 
To hnd typical Nusselt numbers for a simulation we calculate time averaged values for the 
period after secondary saturation. We define the secondary saturation time to be the point 
at which the total kinetic energy of the system reaches a statistically stationary level, and 
we identified it by inspection. To estimate errors we first divide the time average domain 
into 10 bins. We then take the average of each bin and calculate the standard deviation of 
the individual bin averages from the overall average. 

Figure 1^ shows the thermal and compositional Nusselt numbers calculated in both 
cases, for various values of Pr and r. From this figure we also see that the reasonable 
agreement between 2D and 3D simulations discussed in Section persists in all cases 
studied. 

The primary result of this analysis is that, while gravity waves consistently augment 
thermal and compositional transport in non-layered ODDC, the enhanced fluxes are only 
slightly larger than the transport due to thermal and molecular diffusion alone. To be 
precise, the turbulent transport due to gravity waves decreases with increasing Rq^. For 
the simulations we have run with the smallest values of Rq^ (those closest to the layering 
threshold), the turbulent compositional flux is at most twice that of the flux due to diffusion 
alone and the turbulent heat flux is at most ~ 20% of the diffusive flux. However, at larger 
values of Rq^, closer to marginal stability, the turbulent fluxes drop down to ~ 5 — 10% 
of the diffusive fluxes. Also, critically, as Pr and r decrease, the fluxes decrease as well. 
Consequently, the simulations run at the parameter regime most similar to actual giant 
planetary interiors (Pr = r = 0.003) showed the smallest turbulent fluxes compared to 
other simulations with the same value of r, but with larger Pr and r. 

Another result of our analysis is that gravity wave dominated ODDC is responsible 
for the generation of large scale shear. In all the simulations we have run so far, the 
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main effect of the shear has been to modulate the wave-induced transport through strong 
non-linear interactions with the wave held. One might wonder, however, whether the shear 
could become strong enough in some parameter regimes to trigger a shear instability which 
would then dramatically augment turbulent transport. To evaluate the likelihood of this 
happening, we consider the Richardson number, Ri, which is the ratio of the available 
kinetic energy in the shear to the potential energy needed to cause overturn. In the units of 
this paper, we dehne the Richardson number as. 


Ri (z) = 


I 8u |2 
I dz I 


Pr 






(15) 


where N is the buoyancy frequency, dehned dimensionally as N'^ = where p is 

the background density prohle. The terms u and v are the horizontally averaged x and y 
components of velocity, respectively (for 2D simulations u = 0 everywhere, for all time). To 
calculate the typical minimum Richardson number for a simulation we hnd the minimum of 
Ri( 2 ;) for an individual time step, and then take a time average of Rimin(2;) over the period 
after secondary saturation. A plot of the time-averaged minimum Richardson number of 
the available simulations (Figure]^ shows that Rimin increases as r (or equivalently, Rq^) 
increases. This is due to the fact that by dehnition systems with lower Rq^ have a weaker 
stabilizing compositional stratihcation compared to their unstable thermal stratihcation, 
making them more susceptible to overturning. Also, recall from Figure that simulations 
with higher values of Rq^ have a lower percentage of their total kinetic energy is in 
shearing modes. This Richardson number data therefore suggests that if any shear-induced 
instabilities were to present themselves, they would do so at values of Rq^ that are closer 
to the layering threshold. 
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6. Conclusion 


This study marks the conclusion of a series of papers aimed at fully describing 
the thermal and compositional flux properties of ODBC, throughout the entire ODBC 
parameter regime. Rosenblum et ah ( 2011[ ) laid the groundwork for this series by 
conducting a general survey of the OBBC parameter space and identified that OBBC 


either spontaneously form layers, or remains in a non-layered state, which Mirouh et ah 


(2012) later showed to be dominated by large scale gravity waves. They also showed that 


the critical parameter to making predictions about layer formation is the inverse total flux 
ratio defined as, 

Nu^ + y-Jb (Nut - 1) 


7tot = -rRo 


Nu-r 


1 + (Nut’ — 1) 


(16) 


where is the inverse turbulent flux ratio. More precisely, Rosenblum et ah (2011) 
showed that (similar to the conditions that lead to layer formation in fingering convection 


Radko 2003) layers only form when. 


d^tol 


< 0 . 


(17) 


Next, Mirouh et al. (2012) (Paper I) produced a semi-analytical model for y^^) by 
developing an empirically motivated prescription for Nu^ — 1 given by. 


'pr\ 0.25±0.15 

Nut - 1 = (0.75 ± 0.05) ( ^ 

Rn 


T 


(1 -r) , 


r J R- - 1 

where r is the quantity defined in Equation ([^. They also derived an expression for 
that depends only on parameters that can be calculated through a linear analysis of the 
original governing equations in ([^, 

;2 


_ {wfi) _ (Ar R) Xf Ar tR 

Aturb ~ ~ — rin 


(wT) 


''0 


(Ar -|- tR) + Aj Ar -|- R 


(18) 


where Ar and Ai are the real and imaginary parts of the growth rate of the fastest growing 
mode of the primary OBBC instability and I is the horizontal wavenumber of the fastest 
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growing mode (see their appendix for an explanation on how Ar and Aj are calculated, as 
well as analytical approximations in the limit of small Pr and r). 


Importantly, using their model for 7 ^^^, Mirouh et al.l (|2012|) were able to make 


predictions about the range of where layer formation is possible. The function is 
concave-up with a single minimum, so by identifying the value of Rq^ at which reaches 
its minimum (referred to as possible to identify the region of parameter space 

where layers form (Rq^ < R[^) and where they do not (7?o^ > They found that, 

typically 


R 


-1 


tl = 


R-^ 


Pr2 . 


(19) 


Wood et ah (2013) (Paper II) then presented prescriptions for quantifying the thermal 


and compositional fluxes through layered systems. They found that the thermal and 
compositional Nusselt numbers can be modeled as. 


Nu7 


1 = O.lPrsRal 


Nu,, — 1 = 0.03r ^Pr^Ra! 


0.37 

T 


( 20 ) 


The parameter Rar is the thermal Rayleigh number for layered convection, and is dehned 
as a function of the layer height, H, 


, , ,'H\ C(9\Toz 

RMH) = -y = - 

a ) ktv 


TS\H^ 


( 21 ) 


It is not clear a priori what the value of H should be because in simulations of layered 
systems, layers always gradually merge until only a single interface remains. This suggests 
that some other physical mechanism outside the scope of our model of ODBC determines 
layer height. For now, it is left as a free parameter, much like the mixing length in 


mixing-length theory (see also Moore & Garaud 2015) 


By combining the results from Mirouh et ah (2012) and Wood et ah (2013) with the 


work done in this paper, we therefore dehne the effective diffusivities of temperature and 
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chemical composition for the entire ODBC parameter regime, 




efF,^ 


0.03r ^Pr^Ra^^^ + 1) K 


< R 


-1 




B 


eff,T — 


5 


Ro' > RZ' 


O-lPraRa^ + 1 ] Kt , Rq < Rl 

Kt 


-1 


( 22 ) 


, Rq > -Rl 


The implication of onr hndings from this paper for planetary modeling is that 
non-layered ODBC leads to flnxes that are not signihcantly larger than thermal conduction 
or molecular diffusion. The thermal and compositional transport due to gravity waves is 
particularly small when compared to the increase in fluxes caused by the emergence of 
thermo-compositional layers (where turbulent fluxes can be orders of magnitude larger 
the diffusive transport). Consequently, it is not sufficient simply to know if regions in the 
interior of a giant planet are unstable to OBBC. The fact that layered and non-layered 
OBBC lead to very different transport characteristics means that special attention must be 
paid to calculating the threshold R^^, to determine which type of behavior will manifest. 


The dynamics of gravity-wave-dominated OBBC is not expected to be pertinent to 
stars where typical values of R^^ are in the layering regime. However, it may be important 
in giant planets, where higher values of Pr (compared to stars) indicate a smaller range of 
Rq^ that is OBBC unstable, and a lower layering threshold, R^^- 


In particular, the near-zero luminosity of Uranus suggests that thermal transport 


through the planet’s interior is inefficient (Hubbard et ah 1995). Advances in equation 


of state research (Redmer et al. 2011) lends credence to the idea that convection is being 
inhibited by steep compositional gradients. If this is the case, the inefficient, gravity 
wave-dominated, OBBC discussed in this paper may potentially play a role in Uranus’s 
thermal evolution. 


Also, though it is not yet known whether this phenomenon could produce observable 
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signatures, the intermittent growth of shear layers discussed in this work is a potentially 
signihcant feature of non-layered ODBC in the deep interiors of giant planets. In contrast 
to our simulations where there is symmetry between the x and y directions, global rotation 
may provide a preferred direction for shearing motions, which could induce large scale 
azimuthal flows. 

Future studies of ODBC will involve developing more sophisticated models through 
the inclusion of additional physical processes. The natural next step is to include global 
rotation to analyze the effects of Coriolis forces on layer formation and transport properties. 
This may be of particular interest because the gas giant planets in our own solar system are 
rapid rotators, suggesting that giant planets outside our solar system may be as well. 

The authors thank J. Brown, N. Brummell, G. Glatzmaier and J. Fortney for 
enlightening discussions. P. G. and R. M. are funded by NST-AST-0807672 and NSF-AST- 
1211394. The simulations were run on the Hyades cluster, purchased using an NSF MRI 
grant. 
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Fig. 1.— Figure la shows the vertical component of velocity held during the basic insta¬ 
bility growth phase {t = 2508), with red and blue signifying upward and downward motion, 


respectively. Figure lb shows the vertical component of the velocity held at the saturation 


of the primary instability {t = 2868). For each hgure, = 7.87 and Pr = r = 0.03. 
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(a) 



Fig. 2.— The Figures and 2b show, respectively, the velocity held in the y direction at 
times when the system is dominated by gravity waves {t = 20758) and when the system is 
dominated by shear {t = 24718). Here red and blue represent motion in the positive and 
negative y direction. As with the snapshots in Figure= 7.87 and Pr = r = 0.03. 






22 



Fig. 3.— Total kinetic energy vs. time for various families of modes (see main text for 
detail) for a simulation with Pr = r = 0.03 and Rq^ = 7.87. Shown is the early part of a 
simulation where the total kinetic energy is dominated by modes that are predicted to grow 
the fastest according to linear theory. The mode family (6,5,0) is one of the fastest growing 
mode families. 
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Fig. 4.— As in Figure]^ (same parameters), but now focusing on the long-term evolution of 
the system. After the saturation of the primary instability, large scale gravity waves begin 
to grow and energy becomes more concentrated in larger scale modes as time goes on. Large 
scale shearing mode families (mode families (0,0,1) and (0,0,2)) also grow to large amplitude. 
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Fig. 5.— The top figures show the fiuxes of temperature and composition in terms of the 


Nusselt numbers, which are related to the turbulent fiuxes by the equations in (14). Shown 
in the bottom figure is the total kinetic energy in the large-scale gravity wave perturbations 
(modes from the families (l,0,n)) in green, and the kinetic energy in the background shear 
(modes from the families (0,0,m)) in red, as a function of time. The vertical lines mark times 
at which there are bursts of energy in the shear. The simulation parameters used here are 
Pr = r = 0.03, and Rq^ = 7.87. 
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Fig. 6.— Figures 6a and 6b show the thermal and compositional dissipations vs. time for 
the simulation with Pr = r = 0.03 and Rq^ = 7.87. Included are data from a 3D simulation, 


and two 2D simulations with differing domain sizes. 























Fig. 7. 
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Large-scale gravity wave KE (as fraction of totai) 



(a) 


Horizontal KE (as fraotion of total) 



(b) 


shows the energy in gravity wave families of the form (1,0,n) 


Figure 7a 


as a 


percentage of the total energy in the system. Figure [7b| shows the energy in shearing modes 
families of the form (0,0,m) as a percentage of the total energy. We estimate errors according 
to the method described in Section [H 
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Fig. 8.— Figures and ^ show the Nusselt numbers after secondary saturation for the 
available 3D and 2D simulations (in domains of size lOOd^ or lOOd^) as a function of r (which 
is related to Rq^ by Equation]^. The secondary saturation level decreases as r increases 


and as Pr and r decrease. 
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Fig. 9.— The plot shows the time- and horizontally-averaged minimnm Richardson nnmber 
as a fnnction of Rq^- Inclnded are data from all available 2D and 3D simulations (in domains 
of size lOOd^ or lOOd^). 
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